clear all
close all
M=10; %阵元个数
% c=1500; %声速
T=50;
snr=10; % 信噪比
tt=1:T;
L=length(tt); % 快拍数
grid_interval=2; % 网格间距，2<=grid_interval<=4
Azimuth_grid=-90:grid_interval:90; % all hypothetical angles所有假设的角度

timetol=0;iter=0;
TEST=10; %随机试验统计次数
for i=1:TEST
% 角度集：
AL=[ -30 0 20 60 ];
%AL=[ -15.7 -5.3 3.27 11.1 ];
K=size(AL,2); % 信源数量
%% ULA array manifold
A=[];
for m=1:M
for i=1:K
A(m,i)=exp(1j*pi*(m-1)*sin(AL(1,i)*pi/180));
end
end
%% the ULA array completed manifold matrix 
for m=1:M
for ang= 1:length(Azimuth_grid)
Az(ang)=Azimuth_grid(ang)*pi/180;
Aa(m,ang)=exp(1j*pi*(m-1)*sin(Az(ang)));%最后的导向矢量矩阵
end
end
%% 1.Uncorrelated signal : Guassian
S=(randn(K,L)+1j*randn(K,L)); % 产生信号,等功率
%
% Vj=diag(sqrt(( 10.^(snr/10) )/2));
Vj=diag(sqrt( 10^(snr/10)./diag(1/L*(S*S') ) ) );
S=Vj*S; % 产生信号
noise=sqrt(1/2)*(randn(M,L)+1j*randn(M,L)); % 加噪声
Y=A*S+noise;

%% SBL procedure:
% SBL initialize:

[mu,dmu,k,gamma] = sparse_learning(Aa,Y,1,500,1,0,1);


end

 

%% 真实DOA:
plot(Azimuth_grid,gamma)
xlim([-90 90])
xlabel('方位角(度)','fontsize',12);
ylabel('幅度(dB)','fontsize',12);
legend('SBL')



function [mu,dmu,k,gamma] = sparse_learning(Phi,T,lambda,iters,flag1,flag2,flag3)
% *************************************************************************
% 
% *** PURPOSE *** 
% Implements generalized versions of SBL and FOCUSS for learning sparse
% representations from possibly overcomplete dictionaries.
%
%
% *** USAGE ***
% [mu,dmu,k,gamma] = sparse_learning(Phi,T,lambda,iters,flag1,flag2,flag3);
%
%
% *** INPUTS ***
% Phi       = N X M dictionary
% T         = N X L data matrix
% lambda    = scalar trade-off parameter (balances sparsity and data fit)标量权衡参数（平衡稀疏性和数据拟合）
% iters     = maximum number of iterations
%
% flag1     = 0: fast Mackay-based SBL update rules
% flag1     = 1: fast EM-based SBL update rule
% flag1     = 2: traditional (slow but sometimes better) EM-based SBL update rule
% flag1     = [3 p]: FOCUSS algorithm using the p-valued quasi-norm
%
% flag2     = 0: regular initialization (equivalent to min. norm solution)
% flag2     = gamma0: initialize with gamma = gamma0, (M X 1) vector
%
% flag3     = display flag; 1 = show output, 0 = supress output
%
% *** OUTPUTS ***
% mu        = M X L matrix of weight estimates
% dmu       = delta-mu at convergence收敛时的 delta-mu
% k         = number of iterations used结束的迭代次数
% gamma     = M X 1 vector of hyperparameter valuesM X 1 超参数值向量
%
%
% *************************************************************************
% 
% *************************************************************************
    

% *** Control parameters ***控制参数
MIN_GAMMA       = 1e-16;  
MIN_DMU         = 1e-12;  
MAX_ITERS       = iters;%iters
DISPLAY_FLAG    = flag3;     % Set to zero for no runtime screen printouts对于无运行时屏幕打印输出设置为零


% *** Initializations ***初始化
[N M] = size(Phi); 
[N L] = size(T);

if (~flag2)         gamma = ones(M,1);    
else                gamma = flag2;  end;   %初始化超参数向量
 
keep_list = [1:M]';
m = length(keep_list);
mu = zeros(M,L);%初始化权重矩阵
dmu = -1;
k = 0;


% *** Learning loop ***
while (1)%表达式的结果非空并且仅包含非零元素（逻辑值或实数值）时，该表达式为 true
    
    % *** Prune things as hyperparameters go to zero ***当超参数变为零时修剪事物
    if (min(gamma) < MIN_GAMMA )
		index = find(gamma > MIN_GAMMA);%寻找超参数向量大于MIN_GAMMA的索引值
		gamma = gamma(index);
		Phi = Phi(:,index);
		keep_list = keep_list(index);
        m = length(gamma);
     
        if (m == 0)   break;  end;
    end;
    
    
    % *** Compute new weights ***%计算新的权重(w0=Γ0的二分之一次方*(Phi*Γ0的二分之一次方的Moore-Penrose广义逆矩阵)*T)
    G = repmat(sqrt(gamma)',N,1);%G变成块矩阵，N*1的sqrt(gamma)'(共轭转置)N*M(m)Γ0的二分之一次方
    PhiG = Phi.*G;%(N*M矩阵) 
    %求(Phi*Γ0的二分之一次方的Moore-Penrose广义逆矩阵)
    E = pinv(PhiG);%求Moore-Penrose广义逆矩阵
    Xi = G'.*E; 
        
    mu_old = mu;
    mu = Xi*T; 
    
    
    % *** Update hyperparameters ***更新超参数
    gamma_old = gamma;
    mu2_bar = sum(abs(mu).^2,2);%对每一行进行求和
    
    if (flag1(1) == 0)
        % MacKay fixed-point SBLMacKay 定点 SBL
        R_diag = real( (sum(Xi.'.*Phi)).' );
        gamma = mu2_bar./(L*R_diag);  
        
    elseif (flag1(1) == 1)
        % Fast EM SBL
        R_diag = real( (sum(Xi.'.*Phi)).' );
        gamma = sqrt( gamma.*real(mu2_bar./(L*R_diag)) ); 
        
    elseif (flag1(1) == 2)
        % Traditional EM SBL
        PhiGsqr = PhiG.*G;
        Sigma_w_diag = real( gamma - ( sum(Xi.'.*PhiGsqr) ).' );%Xi的转置矩阵
        gamma = mu2_bar/L + Sigma_w_diag;
        
    else
        % FOCUSS
        p = flag1(2);
        gamma = (mu2_bar/L).^(1-p/2);
    end;
    
    
    
    % *** Check stopping conditions, etc. ***检查停止条件，
  	k = k+1;   
    if (DISPLAY_FLAG) disp(['iters: ',num2str(k),'   num coeffs: ',num2str(m), ...
            '   gamma change: ',num2str(max(abs(gamma - gamma_old)))]); end;    
    if (k >= MAX_ITERS) break;  end;
    
	if (size(mu) == size(mu_old))
        dmu = max(max(abs(mu_old - mu)));
        if (dmu < MIN_DMU)  break;  end;
    end;
   
end;



   
end